function out = simulator_threshold(theta2,output,thresholdLevel,numSim)
% We simulate for all maturities
tmp      = output.model;
setup.nn = tmp.nn;
setup.nm = tmp.nm;
nm       = tmp.nm;
setupAll           = setup;
setupAll.matSelect = 1:max(tmp.matSelect);
if output.model.modelNum == 1
    [modelAll,errorMes]= solveQTSM(tmp.alpha,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setupAll);
elseif output.model.modelNum == 2
    [modelAll,errorMes] = solveQTSM_constRsmooth(tmp.alpha,tmp.rhor,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setupAll);
end
modelAll.modelNum = output.model.modelNum;
modelAll.factorSignResOn = 0;
modelAll.rhor     = output.model.rhor;

% Unfold the factor dynamics
nx = size(modelAll.gx,2);
[h0_1,h0_2,hx_1,hx_2,sigma,gama0,gamaz,gamax,sigmaz]   = unfoldTheta2_threshold_macro(theta2,nx);

%% Simulation
rng(1,'twister')
numSim         = numSim + 1000;     %adding a 1000 draws for the burnin
shocks         = randn(nx,numSim);
shocksz        = randn(1,numSim);
xSim           = zeros(nx,numSim);
rLag           = 0;
z              = zeros(numSim,1);
dimy           = size(modelAll.g0,1);
ySim           = zeros(dimy,numSim);
modelAll.macros= zeros(numSim,nm);
for s=2:numSim
    if z(s-1,1) >= thresholdLevel
        xSim(:,s) = h0_1 + hx_1*xSim(:,s-1) + sigma*shocks(:,s);
    elseif z(s-1,1) < thresholdLevel
        xSim(:,s) = h0_2 + hx_2*xSim(:,s-1) + sigma*shocks(:,s);
    end
    z(s,1) = gama0 + gamaz*z(s-1,1) + gamax'*xSim(:,s-1)+ sigmaz*shocksz(1,s);
    modelAll.macros(s,:) = xSim(1:nm,s)';
    modelAll.rLag(s,1)   = rLag;
    ySim(:,s) = modelFit(modelAll,xSim(nm+1:end,s),s); 
    rLag = ySim(1,s);
end

% The regressions
ySim = ySim(:,1001:end)';
z    = z(1001:end,1);
h_index = [1,3,6,9,12,15,18,21,24];
for i=1:length(h_index)
    shortRateOneYear  = ySim(:,modelAll.matSelect==12);
    startIdx          = h_index(1,i);
    endIdx            = max(output.model.matSelect);
    maturities        = modelAll.matSelect(1    ,startIdx:h_index(1,i):endIdx);
    yields            = ySim(1:end,startIdx:h_index(1,i):endIdx);  
    shortRate         = yields(:,1);
   
    % Moments for model evaluations
    tmp = ['h',num2str(h_index(1,i))];
    bootOn = 0;
    out.(tmp) = momentsModelEval(yields,shortRate,shortRateOneYear,z,modelAll.macros(:,2),thresholdLevel,maturities,h_index(1,i),bootOn);
end
out.xSim = xSim(:,1001:end);
out.zSim = z(1001:end,1)';

%% Computing impulse response functions in the two regimes
shockSize = 1;
IRFlength = 60;
numSim    = 5000/100
model     = output.model;
meanx_1   = mean(out.xSim(:,out.zSim>=0),2);
meanx_2   = mean(out.xSim(:,out.zSim< 0),2);
meanz_1   = mean(out.zSim(:,out.zSim>=0),2);
meanz_2   = mean(out.zSim(:,out.zSim< 0),2);
% Adding vecC to model
model.vecC = zeros(length(model.A),model.nx^2);
for k=1:length(model.A)
    model.vecC(k,:) = reshape(model.C(:,:,k),1,model.nx^2);
end
[IRFxhr_1,IRFy_1,IRFx_1] = getGIRFs(model,theta2,meanx_1,meanz_1,numSim,IRFlength,shockSize);
[IRFxhr_2,IRFy_2,IRFx_2] = getGIRFs(model,theta2,meanx_2,meanz_2,numSim,IRFlength,shockSize);
out.GIRFs.IRFxhr_1  = IRFxhr_1;
out.GIRFs.IRFy_1    = IRFy_1;
out.GIRFs.IRFx_1    = IRFx_1;
out.GIRFs.IRFxhr_2  = IRFxhr_2;
out.GIRFs.IRFy_2    = IRFy_2;
out.GIRFs.IRFx_2    = IRFx_2;
out.GIRFs.numSim    = numSim;
out.GIRFs.matSelect = model.matSelect;

end